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COMPUTATION  OF  STATE  VECTOR  AND 
TRANSITION  MATRIX  FOR  TRAM 


I.  INTRODUCTION 

1.1  TRAM  (Trajectory  Reconstruction  Analysis  Method) 

TRAM  will  be  a system  of  subroutines  produced  for  SAMTEC  which,  running 
under  the  control  of  a main  program,  will  use  telemetered  inertial  guidance 
data  and  metric  sensor  data  to  estimate  selected  parameters  of  the  missile, 
the  sensors,  the  earth  and  the  atmosphere.  These  parameters  will  Include, 
for  example,  the  position  and  velocity  of  the  vehicle,  coefficients  of  the 
error  models  for  the  guidance  system  and  for  the  metric  sensors.  Normally, 
parameters  of  the  earth  and  atmosphere  will  not  be  estimated.  It  will  be 
possible,  however,  to  propagate  the  effect  of  errors  of  these,  or  any  of 
the  parameters  not  being  estimated,  on  those  being  estimated.  TRAM  will 
be  useful  for  many  purposes,  such  as  mission  planning  and  special  studies, 
but  it  will  be  used  mainly  for  metric  sensor  performance  analysis  and  will 
subsequently  become  a part  of  the  Metric  Integrated  Processing  System  (MIPS). 

1.2  The  Estimation  in  TRAM 

The  parameters  being  estimated  are  called  the  state  vector.  The  estimation 
is  an  iterative  process  which  takes  place  in  two  steps.  The  first  is  the 
optimal  estimation  of  the  state  vector  and  its  covariance  matrix  at  any 
given  time  using  all  the  data  up  to  and  including  that  time.  The  second 
step  is  to  compute  the  optimal  estimate  of  the  state  vector  and  its  covariance 
matrix  at  any  time  using  all  the  data.  The  first  step  is  called  filtering. 

The  second  is  called  smoothing. 

The  filtering  in  TRAM  will  be  done  by  the  Carlson  square  root  adaptation  of 
the  extended  Kalman  filtering  algorithm  [1].  The  Kalman  algorithm  Is  applic- 
able to  linear  systems  and  in  this  application  the  equations  are  linearized  '{ 
about  a nominal  state  vector. 
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There  will  be  two  options  for  smoothing,  dependent  upon  the  segment  of  the 
trajectory  where  the  data  is  taken  from.  In  powered  flight  a fixed  Interval 
smoother  will  be  used,  and  during  free  fall,  the  smoothing  will  be  done  by 
retrograde  integration.  Figure  1 gives  a functional  flow  of  the  filtering 
and  smoothing  TRAM  will  perform. 

1.3  The  Purpose  of  this  Study 

The  calculations  within  the  dashed  box  in  Figure  1 were  studied  in  this 
effort.  The  nominal  state  vector  Is  computed  by  Integrating  the  position 
and  velocity  from  one  data  time  to  the  next.  With  the  nominal  state  vector, 
the  transition  matrix  can  be  computed.  The  computation  of  the  transition 
matrix  will  be  by  numerical  Integration  also.  Since  the  Intervals  between 
data  times  are  normally  small  (.1  sec  or  less),  and  the  state  vector  and 
transition  matrix  are  needed  at  each  data  point,  computation  times  would  be 
excessive.  Because  of  this,  more  efficient  algorithms  were  required.  These 
were  devised  and  numerical  experiments  were  made  to  test  their  execution 
time  and  accuracy. 

Therefore,  the  main  purpose  of  this  effort  was  to  establish  efficient  algo- 
rithms and  subroutines  that  will  compute  (1)  the  state  vectors  and  transition 
matrices,  (2)  numerical  derivatives  and  (3)  numerical  Integrals  that  are 
needed  for  TRAM. 
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II  CONCLUSIONS 

2.1  The  Three  Areas  of  Study 

Three  problems  were  studied.  One  was  the  size  of  the  integration  step 
which  would  give  adequate  accuracy  when  integrating  the  equations  of  motions, 
the  transition  matrix  and  its  inverse.  Another  was  whether  or  not  numer- 
ical, rather  than  analytic  partials  could  be  used  in  computing  the  transi- 
tion matrix  and  its  inverse.  The  third  was  whether  or  not  Spline  inter- 
polating polynomials  could  be  used  to  obtain  values  of  position,  velocity, 
the  transition  matrix  and  its  inverse  within  the  intervals  of  Integration. 

The  following  results  for  the  free-fall  segment  of  the  trajectory  are  given. 
They  do  not  necessarily  apply  to  powered  flight  or  reentry. 

2.2  The  Conclusions  Concerning  Numerical  Integration 

Numerical  integration  formulas  can  be  divided  into  two  classes,  one  step 
and  multi  step  methods.  The  one  step  methods  use  only  the  Information  at 
the  current  time  step  to  integrate  to  the  next  time  step.  Some  algorithms 
compute  values  of  the  variables  being  integrated  within,  or  at  the  end  of 
the  current  interval  being  integrated  over.  Methods  which  are  of  this  class 
are  the  Euler  method,  the  modified  Euler  method  and  the  Runge-Kutta  methods. 
There  are  many  Runge-Kutta  methods  differing  because  of  the  order  of  their 
accuracy  and  because  of  the  specific  constants  used  in  the  algorithms.  The 
Euler  method,  is  less  accurate  than  the  modified  Euler  which  is,  in  turn,  less 
accurate  than  any  of  the  Runge-Kutta  methods.  However,  the  Euler  method  takes 
less  time  than  the  modified  Euler,  which  In  turn,  takes  less  time  than  any 
of  the  Runge-Kutta  methods. 

The  multi  step  methods  are  numerous.  They  use  Information  at  past  steps  to 
integrate  from  the  current  step  to  the  next.  All  of  these  methods,  however, 
require  that  not  only  the  function  being  integrated  but  its  derivates  up  to 
the  order  of  the  algorithm  be  continuous.  This  is  not  the  case  in  this 
application,  and  therefore,  these  methods  cannot  be  used. 
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Results  from  numerical  experimentation  have  determined  that  the  fourth  order 
Runge-Kutta  gives  the  best  accuracy  and  computation  time  when  compared  with 
the  Euler  and  modified  Euler  methods.  It  has  further  been  determined  that 
a step  size  of  five  (5)  seconds  may  be  used  when  Integrating  the  equations 
of  motion,  the  transition  matrix  and  its  inverse.  These  results  are  valid 
because,  even  though  the  Runge-Kutta  takes  more  time  per  step  than  the 
Euler  or  modified  Euler  methods  Lake,  the  step  size  can  be  appreciably 
larger  because  of  the  much  greater  accuracy  of  the  Runge-Kutta  algorithm. 
Another  reason  for  the  larger  step  size  is  that  during  free-fall,  acceleration 
is  slowly  varying.  During  powered  flight  the  acceleration  is  much  greater,  and 
a smaller  step  size  will  be  needed. 

2.3  Conclusions  Concerning  Numerical  Partial  Di fferentation 

The  differential  equations  of  motion  need  to  be  differentiated  with  respect 
to  position  and  velocity  in  order  to  obtain  the  transition  matrix  and  its 
inverse.  This  study  has  shown  that  these  partial  derivatives  can  be  com- 
puted numerically,  rather  than  analytically,  with  a central  difference 
formula.  The  study  was  made  by  comparing  the  results  of  the  numerical 
differentiation  with  the  analytic  differentiation. 

If  the  analytic  differentiation  is  done  without  mathematical  errors,  and 
the  evaluation  of  the  formulas  is  done  with  infinite  precision,  then  the 
values  would  be  exact.  However,  analytic  partial  differentiation  normally 
becomes  very  detailed  and  numerous  terms  are  developed.  As  a result,  chances 
of  error  are  high  and  very  thorough  checking  is  needed.  Moreover,  if  there 
are  a large  number  of  terms,  round  off  errors  can  degrade  the  accuracy  of 
the  computation. 

Often,  numerical  partial  differentiation  is  simple  to  program  compared  with  the 
analytic  approach.  Once  a subroutine  to  evaluate  the  function  has  been 
developed,  all  that  is  required  to  obtain  the  partial  derivative  with  respect 
to  a given  variable  is  to  evalute  the  function  for  two  different  values  of 
the  variable,  leaving  all  the  other  variables  the  same;  then,  form  a quotient 
of  the  difference  of  the  values  of  the  function  divided  by  the  difference 
of  the  two  values  of  the  variable.  Although  this  is  easy  to  perform  on  a 


computer,  there  are  truncation  and  round  off  errors  associated  with  this 
method  which  the  analytic  approach  does  not  have.  The  reduction  of  these 
errors  to  an  insignificant  value  was  one  of  the  goals  which  this  study  achieved. 

There  is  another  advantage  to  obtaining  partial  derivatives  numerically  and 
that  is  It  is  much  more  flexible.  The  acceleration  calculations  depend  upon 
the  model  of  gravity  acceleration,  of  course.  There  are  many  gravity  models 
and  even  the  form  of  them  can  change.  To  write  down  the  partial  derivatives 
of  all  of  them  would  be  a monumental  task  indeed.  However,  parti als  can 
be  obtained  relatively  simply  by  numerical  differentiation. 

2.4  Conclusions  on  Using  Spline  Interpolating  Polynomials 

The  position,  velocity,  transition  matrix  can  be  integrated  accurately  at 
step  sizes  up  to  five  (5)  seconds  during  free-fall,  but  the  values  of  these 
quantities  will  be  needed  at  every  tenth  second  or  at  smaller  intervals. 

In  this  study  it  was  determined  that  values  at  the  finer  mesh  can  be  obtained 
satisfactorily  by  interpolation  with  Spline  polynomials. 

Spline  polynomials  are  constructed  such  that  the  function  and  a specified 
number  of  its  derivatives  fit  data  at  the  end  points  of  an  interval.  In 
the  case  of  the  state  vector,  the  position,  velocity,  and  acceleration  are 
given  at  the  end  points.  For  each  of  the  three  variables  this  specifies 
six  conditions,  and  therefore,  a quintic  polynomial  is  determined.  For 
the  transition  matrix  and  its  Inverse,  the  function  and  its  derivative  are 
given  at  the  end  points  of  the  Interval.  This  specifies  four  conditions, 
and  therefore,  a cubic  is  determined. 

This  study  has  shown  also  that  the  velocities  in  the  state  vector  can  be 
obtained  satisfactorily  at  the  finer  mesh  by  evaluating  the  derivatives  of 
the  Spline  interpolating  polynomials  that  have  been  derived  for  position. 

2.5  Suitnary  of  Results 

The  amount  of  computation  time  can  be  reduced  appreciably  by  integrating 
the  state  vector,  the  transition  matrix  and  Its  inverse  over  a five  second 
interval  and  then  obtaining  these  quantities  at  the  finer  Intervals  of  the 
data  rate  (0.1  second  or  less)  by  using  Spline  Interpolating  polynomials. 
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It  has  been  determined  that  the  accuracy  of  these  algorithms  Is  sufficient 
for  the  estimation  performed  In  TRAM. 

It  has  been  determined  also  that  the  partial  derivatives  needed  to  calculate 
the  transition  matrix  can  be  computed  by  numerical  partial  differentiation 
with  accuracy  sufficient  for  TRAM.  Computation  of  derivatives  numerically 
will  Improve  computation  times  and  allow  for  greater  flexibility. 

Therefore,  the  purpose  of  this  effort  has  been  accomplished  and  algorithms 
and  subroutines  have  been  developed  that  will  compute,  (1)  the  state  vectors 
and  transition  matrices,  (2)  partial  derivatives,  and  (3)  the  numerical 
integration  needed  for  TRAM. 
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MATHEMATICAL  DESCRIPTION 


3.1  The  Mathematical  Description  of  the  Integration  Algorithm 

The  differential  equations  of  the  motion  of  a vehicle  are  given  by 

p(t)  * v(t) 
v ( t ) = a(p(t),  v( t ) ) 

(3.1.1) 

with  the  Initial  condition  of 
P(o)  = pQ 
v(o)  = vQ 

where  p(t)  is  the  position  at  t,  v(t)  is  the  velocity  at  time  t and 
a(p(t),  v ( t ) ) is  the  acceleration.  The  acceleration  has  the  form 


a(p(t),  v (t ) ) = ag(p(t) ) + ar(p(t),  v(t) ) + ad(v(t)} 

(3.1.2) 


Where  ag  Is  the  gravitational  term  and  Is  a function  of  position,  ar  Is  the 
acceleration  which  would  be  sensed  because  the  equations  are  In  a rotating 
coordinate  system,  and  is  the  drag  and  lift  term.  In  the  case  under 
study  free  fall, 

ad  = °- 


The  differential  equations  of  the  transition  matrix  are  given  by 


Q(t)  = F(t)Q(t) 
Q(0)  = I 


(3.1.3) 


Q Is  a square  matrix  of  order  6 and  F Is  determined  by  differentiating  equations 
(3.1.1)  with  respect  to  position  and  velocity.  F is,  therefore,  equal  to 


F ■ 


3A(j>jV? 
<ipT 


I 

3v 


(3.1.4) 
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where  each  partition  of  F Is  a 3 x 3 matrix. 


Both  (3.1.1)  and  (3.1.3)  are  Integrated  by  the  same  algorithm  and,  therefore, 
the  equations  are  written  in  general  form 


X(t)  = f (X,t) 

X(o)  - XQ 

and  the  algorithm  is  applied  to  the  general  form. 


(3.1.5) 


The  method  used  is  the  standard  fourth  order  Runge-Kutta  method  [2].  It 
is  described  as  follows. 

Let  the  values  of  the  variables  be  given  determined  at  time  t . Therefore, 
Xn  is  known.  The  values  of  the  variables  at  time  tn+^  are  then  determined 
by  these  equations. 

K1  = hf<V  V 

K2  - hf(Xn  + Kr  tn  + h/2) 

K3  = hf(Xn  + 1/2  K2,  tn  + h/2) 
k4  = hf(xn  + k3,  tn  + h) 

Xn+1  = Xn  + 6 ^K1  + 2K2  + 2K3  + K4^ 

In  the  formulas  above  h is  the  step  size  and  the  function  f is  the  right 
hand  side  of  the  differential  equations.  X represents  the  unknown  variables 
being  integrated. 
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3.2  Algorithms  for  Numerical  Partials 

From  equation  (3.1.4)  it  is  seen  that  the  partials  of  acceleration  with 
respect  to  position  and  velocity  are  needed.  Equation  (3.1.2)  gives  the 
form  of  the  acceleration,  with  a^  = 0. 

The  gravity  ar  is  the  rotation  term  and  has  a very  simple  form  which  can 
be  differentiated  analytically  with  ease. 

The  gravity  term  ag  may  vary  markedly  In  Its  form  and  the  number  of  terms 
it  has.  This  leads  to  complicated  analytic  partials  and,  therefore,  the 
numerical  approach  was  used.  The  equation  can  be  derived  rimply  from  a 
Taylor  series  expansion. 
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Let 


XQ  + AX  = X1 


AX 


'-1 


ffX,)  = f(XQ)  + AXf  (XQ)  + f"(X0)+  ^ fM,(n) 


where 


xo  1 xi 


58  f(x0)  - AXf'(x0)  + f"(x0)  - 


Tj-f'U) 


where 


X-1  < 5<  XQ 


Subtracting  the  second  equation  from  the  first,  we  have 

3 

f(X,)  - ftX.,)  » 2AXf'(X0)  + ^y  [f'"(n)  + f"'(0] 
Therefore 

f(x1)  - f(x_1)  ax2  rrM(n)  + f "(e)l 

ZHX " f (Xn)  + -T3T 
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(3.2.2) 


(3.2.3) 


(3.2.4) 
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inis  is  tne  central  difference  formula  tor  tne  derivation  of  fix;  at  xQ. 

The  truncation  error  Is  given  by  the  second  term  on  the  right  hand  side  of 
the  equation.  The  trucatlon  error  Is  seen  to  be  of  order  (AX)2.  Therefore, 
the  formula  used  to  compute  the  partial s Is 

9Aq(X)  Ag(X  + AX)  - Ag(X-  AX)  (3  2 4 

9X  = 2ax 

where  the  value  of  AX  was  determined  to  be  five. 

3.3  Spline  Interpolation 

Spline  interpolating  polynomials  are  used  to  obtain  values  of  the  state 
vector  and  the  transition  matrix  within  the  five  second  interval  of  Inte- 
gration. The  state  vector  uses  a quintic  Spline,  and  the  transition  matrix 
and  its  inverse  uses  a cubic  Spline.  The  development  that  follows  Is 
general  and  therefore,  applies  to  both  cases. 

Suppose  a function  f(t)  and  its  derivatives  up  to  and  including  nth  order 
are  given  over  an  interval  t < t < t^.  Because  of  Improvements  In  the 


round  off 

[0,  1]. 

errors. 

the  interpolating  polynomials  are  developed  over  the  Interval 

Let 

u = 

(t-ta)/T 

where 

T = 

*b  - V 

Define 

g(u)  = 

f(t  + uT)  0 < u < 1 

Q 

Then 

g'(u) 

= T f'(t  + uT) 

Q 

g"(u) 

- T2  f"(ta  + uT) 

a 

g^(u)  = Tn  f^(t  + uT) 

fl 


where  the  prime  denotes  differentiation  with  respect  to  the  variable  In  the 
parenthesis. 
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Let  the  spline  for  g(u)  on  the  Interval  [0,  1]  be  g (u);  and  f (t)  be 
defined  as  follows: 

f(t)  = g((t  - ta)/T) 

Then 

f‘(t)  = }g'((t  - ta)/T) 
f(n)(t)  = ~ g(n)((t  - ta)/T) 

What  is  given  is  f (t  ) , f'(t  ),...,  f^n-1^(t  ),  f(th),  f'(th),..., 

/ n i \ « o au  u 

f'  ;(tb).  From  these  the  following  is  computed. 

g(0)  = f(tfl) 
g*(o)  = Tf ■ (ta) 

g(n_1 ) (0)  = Tn-1f(n“1)(ta) 

g(D  = f(tb) 
g’(D  = Tf*(tb) 

g(n-D(l)  = Tn-lf (n-1 ) 

From  these  it  is  possible  to  get  the  coefficients  of  the  normalized  Spline 
on  [0,  1]: 

g(u)  = Ag  + A^u  + ...  + A2n_i  u ^ 
by  solving  the  system  of  linear  equations 

g(0)  * g(0) 

g’(0)  - g ' (0) 


The  first  three  equations  give  the  values  for  AQ,  A-j , and  The  values 
for  A5,  A4  and  A^  are  determined  by  the  last  three  equations.  This,  then, 
determines  the  Spline  polynomial  for  position  in  the  state  vector.  The 
velocity  is  determined  by  differentiating  this  polynomial  and  dividing 
the  derivative  by  T.  This  would  be  done  for  each  variable,  x,  y and  z. 


The  equations  for  the  transition  matrix  or  its  Inverse  would  be 

% = g(0) 

A}  = g'(o) 

^3  + ^2  * ^1  * ^0  * 9(1) 

3Aj  + 2A2  + A,  = g* (1 ) 

These  equations  can  be  easily  solved  and  the  coefficients  for  the  cubic 
are  thereby  determined.  A cubic  Spline  is  needed  for  each  of  the  36 
elements  of  the  transition  matrix  and  its  inverse. 
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I ■' 


Then  in  general,  f is  the  nth  order  Spline  segment  on  [tfl,  tb]  approximating 
f,  where 


f(t)  = g(u) 


t - t. 


f’(t)  = } g’(u) 


t - t. 


= i ;("-d 


(t)  = ^ g (u) 


t - 1. 


3.4 


The  Generalized  Transition  Matrix 


The  solution  of  the  differential  equation 


X(t)  * F(t)X(t)  + G (t)u(t) 


x(to>  ■ xo 


is  given  by  [3]. 


X (t)  = (Jl  (t,  tQ)  [XQ  + /t  $(t0,  o)G(o)  u(o)dcx] 

to 


where  (J>  (t,  tQ)  satisfies  the  differential  equation 


f (t,  tQ)  - F (t)0(t,  tQ),  <1  (tQ,  t0)  * 


Now 


<>  (t,  t0)p-Uu  tQ)  - I 


Therefore 


i (t*  t0)  V + t ) . o 


and 


<*•  ‘o’  ■ ‘„W‘>  ‘o’  ‘o’ 
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♦"’<*•  t0)  • - ♦"'(t.  t0)F(t) 

♦'’('o’  ‘0>  ’ 1 

Now  let  u(t)  = 0,  and  let  XQ  be  arbitrary. 

Then 

x ( t)  - <J>(t,  t0)X0 
X(s)  = <J)(s.  tQ)X0 

But  also 

X(s)  = $(s,  t)X(t) 

Therefore 

<Ms,  tQ)  = (J)(s,  t)<p(t,  t0) 

Now  let  s = t 

o 

I * ¥ (tQ,  t)(J>(t . tQ) 

and  so 

<P(to.  t)  = f](t,  tQ) 

From  the  equation 

<Ms.  t)  = (p(s,  t0)(p(t0,  t) 

It  is  clear  that 

I 1 

<KS’  t)  = t(s.  t0)(p-1(t,  t0)  (3.4.1 ) 

Since  s and  t are  arbitrary,  (3.4.1)  is  a way  to  calculate  a transition 
matrix  from  one  arbitrary  point  to  another. 
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IV 


DESCRIPTION  OF  SUBROUTINES 


4.1  The  Capabilities  of  the  Subroutines 

With  the  mathematics  set  forth,  it  is  now  possible  to  describe  how  It  was 
implemented.  There  were  two  subroutines  developed.  INTCON  integrates  the 
state  vector  and  the  transition  matrix,  and  then  calculates  the  coefficients 
for  the  interpolating  Spline  polynomials.  It  uses  the  Runge-Kutta  inte- 
gration scheme,  and  when  the  transition  matrix  Is  Integrated,  the  partial 
derivatives  in  the  F matrix,  are  calculated  numerically.  Moreover,  It 
computes  the  Spline  polynomial  coefficients  as  described  In  Section  3.4. 

The  second  subroutine  developed  was  INTERP.  It  calculates  the  Inter- 
polated values  for  the  state  vector  and  the  transition  matrix.  The  tran- 
sition matrix  is  from  one  arbitrary  point  to  another.  The  mathematics 
of  Sections  3.3  and  3.4  have  been  employed  in  the  algorithms  in  this 
subroutine.  INTERP  must  be  called  after  INTCON. 

4.2  Subroutine  INTCON 

The  subroutine  INTCON  calls  two  subroutines  INT  and  CON.  INT  performs  the 
Integration  and  CON  calculates  the  coefficients  of  the  Spline  polynomials. 

4.3  Subroutine  INTERP 

The  subroutine  INTERP  first  calculates  the  Interpolated  values  of  the 
state  vector  for  an  arbitrary  time  t^ . It  then  calculates  the  transition 
from  the  last  time  til  to  the  current  ti . It  does  this  by  finding  the 
inverse  of  the  transition  matrix  from  the  beginning  of  the  Interval  to  t^ 

Then  it  finds  by  interpolation  the  transition  matrix  from  the  beginning 
of  the  interval  to  t^  and  multiples  the  two  matrices  together. 

4.4  Other  Subroutines 

These  two  subroutines  call  others.  Altogether  there  were  fourteen  subroutines 
developed. 
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V. 


NUMERICAL  VERIFICATION 


5.1  Verification  of  the  Numerical  Integration 

The  first  steps  taken  to  determine  which  of  the  one  step  methods  would  be 
adequate  were  to  run  the  three  methods,  the  Euler  method,  the  modified  Euler 
method  and  the  Runge-Kutta  method,  over  a large  segment  of  the  trajectory 
and  see  if  adequate  answers  could  be  obtained.  Only  the  Runge-Kutta  gave 
answers  which  were  accurate  enough. 

With  the  method  of  integration  selected  it  was  then  possible  to  determine 
an  optimum  step  size  which  would  give  the  best  accuracy.  A step  size  of 
one  second  and  a step  size  of  five  seconds  gave  answers  which  agreed  to  more 
than  0.001  foot  in  position  and  to  0.00001  foot  per  second  in  velocity  when 
integrated  over  an  interval  of  1500  seconds.  Consequently,  the  five  second 
step  size  was  chosen. 

The  magnitudes  of  the  round  off  and  transition  errors  caused  by  the  numerical 
integration  are  well  below  errors  from  other  sources  and  therefore,  will  not 
affect  the  estimation  capability  in  TRAM.  The  total  CPU  time  to  Integrate  the 
state  vector,  the  transition  matrix  and  its  Inverse  1500  seconds  was  two 
minutes  and  32  seconds.  This  is  an  acceptable  execution  time. 

5.2  Verification  of  the  Accuracy  of  the  Numerical  Partials 

The  transition  matrix  is  Initialized  to  the  identity  matrix  every  five  seconds 
and  then  integrated  for  a five  second  span.  To  bound  the  error'  caused  by  the 
use  of  the  numerical  partials,  this  was  done  twice,  once  using  analytic  partiais 
and  once  using  numerical  partials,  and  the  difference  between  the  two  transition 
matrices  was  computed.  This  difference  matrix  is  then  multiplied  by  a vector 
which  has  as  its  elements  the  largest  error  which  can  occur  at  any  step  In 
position  and  velocity.  The  product  of  the  matrix  times  the  vector  gives 
the  maximum  error  that  can  occur  at  each  step  caused  by  the  use  of  numerical 
partials.  If  the  resulting  vector  is  multipled  by  the  number  of  steps,  a 
bound  for  the  total  error  due  to  numerical  partials  Is  obtained. 
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In  the  error  matrix  which  was  obtained  by  subtracting  the  transition  matrix 

computed  using  analytic  partials  from  the  transition  matrix  using  numerical 

-14 

partials  no  term  was  larger  than  10  . This  meant  that  the  terms  which  propa- 

gates position  and  velocity  errors  into  position  and  the  terms  which  propa- 

-14 

gate  position  and  velocity  errors  into  velocity  were  all  less  than  10 
Further,  the  maximum  TRAM  estimation  error  which  can  be  made  In  position  at 
any  step  is  less  than  100  feet  and  in  velocity  It  is  one  foot  per  second. 

Since  the  data  could  be  taken  at  one  tenth  second  for  1500  seconds,  their 
would  be  15000  steps.  These  numbers  provide  the  following  bounds  for  the 
errors  caused  by  the  use  of  numerical  partials. 


< 

15000 

X 

1 

O 

lO"14^ 

X 

( 102  1 

£iv  I 

^ 10'14 

10-14  ) 

w / 

e-j  p < 15.15  x 10"9  feet 

elv  £ 15.15  x 10"9  feet/sec 


where  e^p  is  the  error  in  position  caused  by  the  use  of  numerical  partials 
and  e^v  In  the  error  In  velocity  caused  by  the  use  of  numerical  partials. 
Clearly,  these  errors  are  insignificant. 

5.3  Verification  of  Accuracy  of  Spline  Interpolation 

To  determine  what  the  accuracy  of  the  Spline  Interpolation  was,  the  trajectory 
was  integrated  to  1500  seconds.  At  the  beginning,  at  mid  trajectory  and  at 
the  end  of  the  trajectory.  Spline  polynomials  were  fit  over  the  five  second 
interval.  The  state  vector,  the  transition  matrix  and  Its  Inverse  were 
approximated  when  analytic  and  also  when  numerical  partials  were  used. 

The  error  In  position  in  the  state  vector  was  always  less  than  10”4  feet  and 
the  error  In  velocity  was  less  than  10"5  feet  per  second.  In  the  transition 
matrix,  the  term  which  propagates  errors  In  position  Into  position  Is  less 

_g 

than  10  , and  the  term  which  propagates  velocity  error  Into  position  Is 

less  than  10'7.  The  term  which  propagates  position  error  Into  velocity 
was  less  than  10”^,  and  the  term  which  propagates  velocity  error  Into 
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velocity  was  10  . Since  position  estimation  errors  are  less  than  100  feet 
per  step,  and  velocity  estimation  errors  are  less  than  one  foot  per  step  and 
there  are  15000  steps,  the  error  terms  are  bounded  by 


1 10'9 

10'7^ 

( >°2 

\£2  •>) 

< 

^ 10' 11 

'°~9I 

X 

l’  i 

x 15000 

where  e£p  and  egv  are  errors  In  position  and  velocity  respectively  which  are 
caused  because  of  Spline  interpolation. 


Therefore 


£2p 
e2  v 


<_  3 x 10~3  feet 

< 3 x 10'5  feet/sec. 


5.4  Numerical  Error  Summary 


For  the  state  vector,  the  two  sources  of  error,  one  from  numerical  Inte- 
gration and  the  other  from  Spline  interpolation  cause  the  following  total 
error 


it 

< 

10'3  ♦ 

1 

O 

= .0011  feet 

P 

nu 

< 

-5 

10  3 ♦ 

10'5 

= 2 x 10"^  feet/second 

where  np  and  ny  are  position  and  velocity  errors  because  of  ..umerical 
approximation.  Both  of  these  errors  are  insignificant. 


For  the  transition  matrix  the  following  total  error  bounds  are  valid. 

Ept-elp  + c2p  = 2 x 10’8  + 3 x 10'3  < 3.1  x 10'3 


Evt  - elv  + e2v  = 2 x 10"8  + 3 x 10*5  < 3,1  x 10"5 
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Both  of  these  errors  are  insignificant.  It  may  be  concluded,  therefore, 
that  the  errors  due  to  the  numerical  approximations  discussed  here  will  in 
no  way  impair  the  estimation  accuracy  in  TRAM  and  further,  it  may  be  con- 
cluded that  these  algorithms  will  greatly  improve  the  speed  and  flexibility 
of  the  computation  of  the  estimation. 


VI 


DOCUMENTATION  OF  SUBROUTINES 


6.1  GENERAL  COMMENTS  ON  DOCUMENTATION 

The  documentation  which  follows  is  intended  to  be  sufficient  enough  so  that 
the  capabilities  of  these  subroutines  can  be  understood,  used  and  modified. 

The  function  that  the  subroutine  performs  is  stated  first.  This  is  followed 
by  a description  of  the  input  and  output  variables.  The  restrictions,  if  any, 
under  which  the  subroutine  must  be  used  come  next,  and  then  a list  of  the 
subroutines  this  subroutine  calls  is  given.  The  documentation  also  includes 
a listing  of  the  subroutine  and  a flow  chart. 

The  main  subroutines  are  INTCON  and  INTERP,  but  these  call  others  which  are 
documented  also. 
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FUNCTION: 


INPUT: 


/VAR/ 


OUTPUT: 


/VAR/ 


PH  I NO 


PHIDO 


INTCON  SUBROUTINE  DESCRIPTION 

Integrates  the  state  yector,  the  transition  matrix  and  Its  Inverse 
over  the  interval  of  Integration  and  then  computes  the  coefficients 
for  the  interpolating  Spline  polynomials.  A quintlc  Spline  is 
determined  for  the  state  vector  and  a cubic  Spline  is  determined 
for  the  transition  matrix  and  its  Inverse. 


Nine  dimensional  state  vector  containing  position,  velocity  and 
acceleration.  Only  position  and  velocity  are  needed  as  Input  to 
INTCON. 

Time  at  the  beginning  of  the  Interval  of  integration. 

Start  time  of  the  trajectory 
Step  size  of  the  integration 

Flag  to  Indicate  whether  or  not  the  inverse  of  the  transition 
matrix  should  be  integrated. 


Nine  dimentional  state  vector  which  contains  position,  velocity 
and  acceleration  integrated  to  the  end  of  the  Interval 

6x6  dimensional  array  containing  elements  of  the  transition 
matrix  at  the  end  of  the  interval  of  integration 

6x6  dimensional  array  containing  the  derivatives  of  PHI 

6x6  dimensional  array  containing  the  elements  of  the  inverse 
of  PHI 

6x6  dimensional  array  containing  the  derivatives  of  PHIN 

6x6  dimensional  array  containing  the  elements  of  the  transition 
matrix  at  the  beginning  of  the  Interval  of  Integration 

6x6  dimensional  array  containing  the  derivatives  of  the  elements 
of  PHIO 


PHINO  6x6  dimensional  array  containing  the  Inverse  of  PHIO 

PHINDO  6x6  dimensional  array  containing  the  derivatives  of  PHINO 


XO  9 dimensional  array  containing  the  values  of  the  state  vector  at 
the  beginning  of  the  interval 

XC  6x3  array  containing  the  coefficients  for  the  quintic  Spline 

for  the  state  vector.  The  first  index  the  coefficients  and  the 
second  indicates  the  variable. 

PHIC  4x6x6  array  containing  the  coefficients  for  the  cubic  Spline 

for  the  transition  matrix.  The  first  index  indicates  the  coef- 
ficient and  the  last  two  indicate  the  element  in  the  transition 
matrix 

PHINC  4x6x6  array  containing  the  coefficients  for  the  inverse  of 
PHI.  The  Indices  are  the  same  as  for  PHIC. 

SUBROUTINES  USED:  FMATRX , I NT,  CON 
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INTCON  FLOWCHART 


IF  T - TS 
INITIALIZE 

DIFFERENTIAL  EQUATIONS 


SET  UP  VALUES  FOR 

STATE  VECTOR  AT 
BEGINNING  OF  INTERVAL 


IF  ID  - I 

INITIALIZE 
PHIN,  PHINO, 
PHINDO,  PHIND 


CALCULATE 
PHID  & PHIDO 


INITIALIZE 
F MATRIX 


CALL  I NT 

INTEGRATE  STATE  VECTOR, 
TRANSITION  MATRIX  & 
INVERSE  FORWARD  i 


LEVEL  21.7  ( jan  73  ) 


OS/360  FORTRAN  H 


UAH 


COMPILER  OPTIONS  - NAMfc  = MAIN.OPT=01  .LlNECNT*e>0.SIZE=0000K, 

SOUkCE .EbCDIC.NOL 1ST .NOOECK .LOAD ,NO«AP .NOEOI T • 10. NOXREF 


ISN 

0002 

SUBROUTINE  INI CON 

ISN 

0003 

IMPLICIT  01 al*6(A-H,0-Z) 

ISN 

000A 

Common/ </ AO/ X (9)  .Pm  1(6. 6)  . PM l lit 6, 6)  .PM IN (6. 6)  .PH] NO  1 6. 6)  , 

I PH  10(6.6)  *PRl DO (6.6) «PHlN()<6»6) .PhIn00(6.6)  <X0(9).XC(6<3>. 

? PhIC IA.6.6)  iPmInC IA*6i6|  ■ ISiT.h.ID 

ISN 

ooos 

COmmON/PHIm/F (6.6) 

ISN 

0006 

DIMENSION  OX (6) 

ISN 

0007 

IF (T .NE.TS)  GO  10  5 

ISN 

0009 

CALL  01FFUNI6.T.X.DX) 

ISN 

0010 

DO  3 1*1.3 

ISN 

0011 

3 

X(I*6)*0X(I»3) 

ISN 

0012 

5 

CONTINUE 

ISN 

0013 

UO  10  1*1.9 

ISN 

001A 

10 

XO(I>=X(l) 

ISN 

001S 

CALL  FMATRXU) 

ISN 

0016 

DO  20  1=1.6 

ISN 

0017 

DO  20  J* 1.6 

ISN 

0016 

PmJ  < I , J)*0D0 

ISN 

0019 

1FU.EQ.J)  PM  I ( I . J)  * 1 DO 

ISN 

0021 

Ph 10(1 . J>  *Phl ( 1 . J) 

ISN 

0022 

20 

CONTINUE 

ISN 

0023 

IF ( IO.EO.O)  GO  TO  19 

ISN 

0025 

DO  22  1=1.6 

ISN 

0026 

00  22  J=1 .6 

ISN 

0027 

PHINI r.J> =PMI (1,J) 

ISN 

0028 

PHlNOII.J)=PMl (I,J) 

ISN 

0029 

22 

CONTINUE 

ISN 

0030 

DO  18  J=1 .6 

ISN 

0031 

00  18  1*1,6 

ISN 

0032 

PHINDOU  ,J)  =000 

ISN 

0033 

DO  17  K=l .6 

ISN 

003A 

1 7 

PNINDOI I . J) = -PH I NO (I»KI*F(K*J) .PHlNOO ( I . J) 

ISN 

003S 

18 

PH I NO  1 1 • Jl =PR I NOO ( 1 . J ) 

ISN 

0036 

19 

CONTINUE 

ISN 

0037 

C Al  l DPMI  (PHI  .PH  10) 

ISN 

0036 

DO  2S  1=1.6 

ISN 

0039 

00  25  J=l .6 

ISN 

00A0 

25 

PM  I DO ( I . J) =PmIO< I , J) 

ISN 

OOA  1 

CALL  INT ( X ,DX .PhI .PhIO.PmIN.PHINO. T .M. 10) 

ISN 

00A2 

CALL  CON 

ISN 

00A3 

RETURN 

ISN 

OOAA 

ENO 

•OPTIONS  IN 

EFFtCT* 

NAME*  MAIN. OPT  = 01  .L  lNtCNT=60.SI/t  = O()00K, 

•OPTIONS  IN  EFFtCT*  SOURCE .EBCDIC iNOL 1ST . NODE CK .LOAD.NOMAp.NnEOIT. 10 .NOXREF 


•STATISTICS*  SOURCE  STATEMENTS  * A3  .PROGRAM  SIZE  * 1132 

•STATISTICS*  NO  DIAGNOSTICS  GENERATED 


67K  BYTES  OF  CORE  NOT  USED 


best  available  cm 


Eno  of  compilation 


6.3 


INTERP  SUBROUTINE  DESCRIPTION 


FUNCTION:  Interpolates  the  state  vector  to  time  TI  using  Spline  polynomials, 
obtains  a transition  matrix  form  TI-HS  to  TI  by  interpolation 


INPUT: 


/VAR/ 


/TST/ 

S 

OUTPUT: 


/TST/ 


Interpolation  time 


HS  Difference  between  last  interpolation  time  and  TI 


Time  at  the  beginning  of  the  interval 


ID  Flag  to  indicate  how  the  inverse  of  the  transition  matrix 
is  to  be  computed 

PHIO  6x6  array  containing  transition  matrix  at  the  beginning  of 
the  interval 


6x6  array  containing  elements  of  transition  matrix  from  the 
beginning  of  the  interval  to  time  of  last  Interpolation  TI-HS 


6 dimensional  array  containing  state  vector  Interpolated  to 
time  TI 

6x6  array  containing  the  transition  matrix  from  time  TI-HS 
to  time  TI 


QI  6x6  array  containing  inverse  of  transition  matrix  from  beginning 
of  interval  to  TI-HS 

S 6x6  array  containing  transition  matrix  to  time  TI 
SUBROUTINES  USED:  PVXTRP,  DMINV,  INXTRP,  FIXTRP,  DGMPRD 


' BESFA\ 

Ml}  art 

cm 

ULJi  r\\ 

r rMLHULi. 

Lor  i 

s 

* LEVEL  21.7  ( JAN  73  ) 

OS/360  FORTRAN  M 

DATE 

COMPILES  OPTIONS  - NAME  = MAIN,OPT*01 .LlNEtNT*60»S!ZE*nOOOK. 

SOURCE. EBCDIC. NOLI  ST. NOOECK, LOAD, NO*' AP.NOCOIT, 10. NOXREf 
ISN  0002  SUBROUTINE  I NTERP ( T I . T .0 . HS) 

C SUBROUTINE  INTERPOLATES  POSITION  AND  VELOCITY  TO  TIME  TI  USING  THE 
C SPLINE  POLYNOMIAL  ANO  STORES  IT  IN  Y.  IT  ALSO  COMPUTES  THE  TRANSITION 
C MATRIX  FROM  TI-mS  TO  TI  USING  THE  SPLINE  POLYNOMIALS  ANO  STORES  IT  IN 
C 0. 


ISN 

0003 

IMPLICIT  REAL*8(A-M,0-Z» 

ISN 

0006 

COMMON/VAR/X (V) .PHI (6.6) .PHI0( 6,6) ,PHIN (6.6) .RHINO (6.6) ♦ 

1 PHI  0(6,6) ,PnIU0<6.6> .PHI  NO (6,6) ,PHIND0(6.6) • XO ( 9) , XC (6, 3) . 

2 PHIC(6,6,6) ,PHl NC (6,6.6) .TS.T.H.ID 

ISN 

OOOS 

COMMON  /TST/  S<6,6) .QI (6,6) 

ISN 

0006 

DIMENSION  Y(6), 0(6. 6), LI (6), Ml (6) 

ISN 

0007 

CALL  PVXTRP(Tl.Y) 

ISN 

0008 

TU®  T*hS* 10-6 

ISN 

0009 

IF ( T I .LE . 2D- 1 ) WRITE (6,100)  TI.TU.T 

ISN 

0011 

100 

FORMAT  <2X, • TI ,TU. f .3E20.6/) 

ISN 

0012 

IF(TI.GT.TU)  GO  TO  10 

ISN 

0016 

DO  5 1*1.6 

ISN 

0015 

00  5 J*1 ,6 

ISN 

0016 

0(I.J)*PHlOII*J) 

ISN 

0017 

5 

S ( I ,J) =Q( I , J) 

ISN 

ooia 

IF (TI.EO.T)  RETURN 

ISN 

0020 

10 

CONTINUE 

ISN 

0021 

IFtlO.EG.l)  GO  TO  15 

ISN 

0023 

DO  12  1=1.6 

ISN 

0026 

DO  12  J=1 .6 

ISN 

0025 

12 

om.j)*sii*j) 

ISN 

0026 

N*6 

ISN 

0027 

CALL  DMlNV(OI.N.O.Ll.Ml) 

ISN 

0026 

GO  TO  20 

ISN 

0029 

IS 

TE«TI-HS 

ISN 

0030 

CALL  INXTRP(TE.OI) 

ISN 

0031 

20 

CALL  FIXTRP(TI.S) 

ISN 

0032 

1FITI.LE.2D-1)  WRITE (6.110)  S 

ISN 

0036 

no 

F ORMAT (2X«,S,/,6(2X,6E20,10/)/) 

ISN 

0035 

CALL  DGMRRD (S.OI ,0.6.6. 6) 

ISN 

0036 

RETURN 

ISN 

0037 

ENO 

• OPTIONS  IN  EFFECT*  NAME » MA In.0PT*0 I ,L I NECNT*60 . S I ZE*OOOOK , 

•OPTIONS  IN  EFFECT*  SOURCE .EBCDIC .NOL 1ST .NOOECK .LOAD .NOMAP .NOEDI T . ID. NOXREF 

•STATISTICS*  SOURCE  STATEMENTS  * 36  .PROGRAM  SIZE  • 10<»6 

•STATISTICS*  NO  DIAGNOSTICS  GENtRATEO 


ENO  OF  COMPILATION 


75K  BYTES  OF  CORE  NOT  USED 


6.4 


INT  SUBROUTINE  DESCRIPTION 


FUNCTION:  Integrates  state  vector,  transition  matrix  and  Its  Inverse  from 
T to  T+H 


INPUT: 


OUTPUT: 


9 dimension  array  state  vector  containing  position,  velocity  and 
acceleration 

6x6  array  containing  transition  matrix  at  time  T 

6x6  array  containing  inverse  of  PHI  at  time  T 

Flag  to  indicate  if  inverse  of  transition  matrix  should  be 
integrated 

Step  size  of  integration 


PS  6x4  array  contains  four  intermediate  values  of  the  state  vector 
used  in  the  Runge-Kutta  integration  algorithm  for  integrating 
the  transition  matrix  and  its  inverse 

X 9 dimensional  array  state  vector  containing  position,  velocity 
and  acceleration  at  T+H 

PHI  6x6  array  containing  transition  matrix  at  T+H 

PHIN  6x6  array  containing  inverse  of  PHI  at  T+H 

DX  6 dimensional  array  containing  derivative  of  X at  T+H 
PHID  6x6  array  containing  derivative  of  PHI  at  T+H 
SUBROUTINES  USED:  DIFFUN,  INTQ,  DPHI,  DPHINV 


LtVt'L  21.7  ( JAN  7 3 ) 


OS/360  FORTRAN  H 


DATE 


ISM 


0002 


COMPILER  OPTIONS  - NAM)  = MAlN,OPT=0l.LlNECMT=S0.SI2E=00n0K, 

SOURCE  .EBCDIC  .NOL  1ST  . NOOt C * . LOAD . NO” AP , nOE  D I T , 1 0 , NOXREF 
SUHPOUT  I NT  liT  ( a.DX.PhI  . Pr.  1 G.  PM  1 N.  PHI  NO.  T .H,  I T)  > 

c sum-routine  int  (.ives  the  values  of  position  ano  velocity,  the 

C TRANSITION  MAT”!.  ANO  ITS  utPIVATlVt  AMU,  AS  AN  OPTION*  Trif  INVEPSf 

C of  the  Transition  and  its  derivative  at  timc  t*h.  io  is  a ei.ac  o< 

C ISDlLATtS  *NtT"f<  THE  INVEPSf  G*  Tt*E  TRANSITION  maTpIa  ANO  ITS 
C OtOIVATt  IS  To  oE  COMPUTED.  The  POSITION  ANO  VE . oC l T Y APE  STOfi  0 * ' 
c TrtF  TRANSITION  MATRIX  AND  IIS  lf.Vtf.SE  APt  STOPEP  IN  Pul  ANl>  P*-'0 

c tmf  inverse  of  transition  matpia  and  its  derivative,  ape  stored  in 

C P M J N AND  PMIf'D. 


* . 


IS* 

0003 

IMPLICIT  REAl«8(A-m.O-2< 

ISM 

0004 

DIMENSION  X U> .OX  < l > .PHI  16.6) .PflD 16,6) .PhIn 16.6) ,PMlN0'6. 
I.  P V ( 6 ) ,PS(h.4),Y(6).0I6).0Y(6). 
f.  C ( 3)  .014) 

1SN 

OOOS 

EXTERNAL  DPrtl.OPHlNV 

ISN 

0006 

DATA  C/2*0 . 500 .100  /,  0/ 1 Dl) . 2*200 . 1 00/ 

ISN 

0007 

DO  S 1=1,6 

ISN 

0008 

PV ( I ) =X  < 1 ) 

ISN 

GOON 

5 

PS ( I . 1 ) =X ( 1 ) 

ISN 

0010 

N=6 

ISN 

0011 

DO  TO  1=1.4 

ISN 

U0I2 

IF ( l .EQ. 1 ) GO  TO  IS 

ISN 

001<* 

K = J-1 

ISN 

001S 

DO  VO  J = i .6 

ISN 

0016 

PV  U) =PS 1 J. 1 ) »C (K ) *0  < J) 

ISN 

0017 

10 

PS(J.I)=PV(J) 

ISN 

C 0 1 1 

IS 

CONTINUE 

ISN 

00  1 P 

CALL  OIFFUNIN.T.PV.OYI 

ISN 

0020 

CO  20  J=1 .6 

ISN 

0021 

20 

0(J)=H*0Y<J) 

ISN 

0022 

00  25  J=l .6 

ISN 

0023 

2S 

X ( J) = x ( Jl *0( J) *0( l 1 /600 

Jsn 

0024 

30 

CONTINUE 

ISO 

0025 

CA|  l PIFFuN(N.T.X.DX) 

ISN 

0026 

DO  31  1=1.3 

ISN 

002» 

31 

x ( 1 »6)  =0*.  ( 1*3) 

ISN 

002B 

CALI  INTO  (PS.PHl  .PHltl.T  .H.X.OPHl  I 

ISN 

002S* 

If  <10. tO. 1)  CALL  INTQIPS.PHlN.PfilNO.r.rt.X.DPHlNV) 

ISN 

0031 

»t  turn 

ISN 

0032 

ENO 

♦OPT  TONS  IN  Ff  E t C T • 


NA«t=  MA  lN,OPT=Ql  .L  INELNT=6C.NI /t  -0000*.  . 


•OPTIONS  IN  EFFECT • 


SOURCE .triCOIC .NOLI  ST .NOOECf  .LOAD , NOMAP .N-EL I T , I0.NGXREF 


♦STATISTICS*  S0OPCE  STATEMENTS  = 31  * PROGRAM  S12t  - 1464 

•STATISTICS*  NO  DIAGNOSTICS  generateo 

••••••  END  of  COMPILATION  ••••*•  7SK  bytes  of  CORE  not  useo 


best  available  copy 


31 


6.5 


PVXTRP  SUBROUTINE  DESCRIPTION 


I 


FUNCTION:  Interpolates  state  vector  to  time  TE 
INPUT: 

TE  Time  at  which  the  interpolating  polynomial  is  to  be  evaluated. 
/VAR/ 

XC  6x3  array  containing  the  coefficients  for  the  Spline  interpolating 
polynomial  for  the  state  vector.  The  first  indicates  the  coef- 
ficient and  the  second  indicates  the  variable  being  interpolated. 

OUTPUT: 

Y 6 dimentional  array  containing  the  interpolated  values  of  position 
of  the  state  vector 

ENTRY  FIXTRP 

FUNCTION:  Interpolates  values  of  transition  matrix  to  time  TE 
INPUT: 

TE  Time  at  which  Interpolating  polynomials  are  to  be  evaluated 
/VAR/ 

PHIC  4x6x6  array  containing  coefficients  of  Spline  polynomials.  The 
first  index  Indicates  the  coefficient  and  the  last  two  indicate 
the  variable 

OUTPUT: 

Q 6x6  array  containing  elements  of  transition  matrix  at  time  TE 
ENTRY  INXTRP 

FUNCTION:  Interpolates  inverse  of  transition  matrix  to  time  TE 
INPUT: 

TE  Time  at  which  the  interpolating  polynomial  is  to  evaluated 


/VAR/ 

PHINC  4x6x6  array  containing  coefficients  of  the  Spline  polynomials  for 
the  Inverse  of  the  transition  matrix.  The  first  Index  Indicates 
the  array  and  the  last  two  indicates  the  element. 

OUTPUT: 

S 6x6  array  containing  the  elements  of  the  inverse  of  the  transition 
matrix  for  time  TE 


EVALUATE  6X6  SPLINES 
FOR  VALUES  OF 
TRANSITION  MATRIX 


LEVEL  21.7  « JAN  73  ) 


OS/3bO  FORTRAN  m 


DATE 


ISN  000? 


ISN  uOUJ 

i sn  ooo* 


ISN  «005 
ISN  0006 
ISN  0007 
ISN  O00O 
ISN  000* 
ISN  0010 

isn  con 

ISN  0012 

ISN  0013 
ISN  001** 
ISN  001S 
ISN  0016 


COMPILER  OPTIONS  - NA«t=  MMN.OPnoi  .LINECnT=6O,SI2E=OOO0k. 

SOUwCE . EoCO I C . NOL 1ST .NOOECK .LO AO .NO  - AP .NOE  U I T . ID .NOXPEF 
subroutine  Pvxr«piTt*vi 

C S'JPpOijllNf  PVXTWP(TE.Y)  F XTWAMOi  Alt  O POSITION  A-it,  VELOCITY  TO  TIME  TE 
C US  I M>  A rjMINTIC  SPLlNf.  fOP  POSITION  AND  ITS  DERIVATIVE  F OP  VELOCITY. 

C The  COEFFICIENTS  FOP  Tnf  QUINT!1:  SPUNE  AWE  COMPUTED  IN  SUMwOUT  1 NE 
C INTCON.  INTERPOLATED  POSITION  /.NO  VELOCITY  A«t  SToRt.0  IN  Y. 

implicit  peal»m<a-m.o-2) 

COMMON/ V A P/X (9)  .pm  I (6.6)  .PHlO (6,6)  .PH  I, VI  6, 6 I . PH  I NU (6, 6) . 

1 PhI0«6.6> .Pm  I DO (6.0) , Pul  NO (6. 6) .PnlNOO<6.6) . XO ( * ) . XC 1 1 . 3)  . 
c ^pIC  ( A.o  ifc)  »Ph1  NC  lw.6.6)  « T S • T >h.  It) 

DIMtNSION  Q(6.6) 

DIMENSION  S ( 6 • 6 ) 

DIMENSION  Y ( 6) 

U=(TF-T)/M 
DO  10  1*1.3 
J=  1 ♦ 3 

Y ( I ) = (JAU^UotJAU^XC  ( 1 . I ) ♦ U^U'U'U^XC  <2. 1 > »U«U»U*XC  < 3 . 1 . »U*U*  xC  «<*.!> 
1 »U»XC  «*>.  I ) •XC(6,I ) 

Y ( J)  =5Q0»U*U«U»U»XC  ( 1 . I ) ♦ AO0»U»U»U«XC  <2.  I ) ♦ 3D0»U»U»XC  < j.  I ) 
l ♦dUU»U*xC(<‘.l)«XC(S.I> 

Y ( J) = Y ( J) /H 
10  CONTINUE 
PE TOWN 

EnTPY  F1XThP(TF.OI 

C SJHWOUTINF  FI>TpP(TE.O)  EXTRAPOLATES  PHI  TO  TIMf  TE  USING  A CUBIC 
C SALINE,  the  COEFFICIENTS  FOW  ThE  CUUIC  SPLlNL  AWE  COMPUTED  IN 
C S'JMPOUTlNf  I-TCON.  INTEPPOLATEU  VALUES  OF  PrtI  ARE  STORED  IN  O. 


ISN 

0017 

U= ( TF-T) /H 

ISN 

ooia 

DO  15  1=1.6 

ISN 

001* 

uO  15  J=  1 *6 

ISN 

0020 

15 

9 ( I . J»  =1)«U»U*PHIC  < 1 . I . Jl  ♦U»U»PHIC (2,  I . J)  *U»PH  [ C ( 3 . 1 . J ) »PHlC(A.  1 . J) 

ISN 

0021 

RETURN 

ISN 

0022 

Entry  inxtppite.S) 

C SUUmOu 1 1 NE  I NX  r-w  I Tc. . S ) FXTWAPOLATES  THE  lNVEWSE  OF  PHI  TO  TIME  TE 
C USING  A CUBIC  SPLINE.  THE  COEFFICIENTS  FOW  THL  CUbIC  SPLINE  AWE 
C CO-PUTtD  IN  SUfcwQuTlNF  INTCON.  INTERPOLATED  VALUES  OF  Pm!  INVERSE  APF 
C STORE  0 IN  S. 

ISN  0023  U=(TE-T)/m 

ISN  002A  DO  20  1=1.6 

ISN  0025  DO  20  J=1 .6 

ISN  0026  20  S ( I . J)=l)»l)»u»PHlMC  ( 1 . 1 . J)  ♦U«U*PHlNC  (2. 1 . Jl  ♦I/APHINC  (3. 1 • J)  ♦ 

2 pMlNCCA. I . J) 

ISN  V02  7 RETURN 

ISN  002t*  END 


•OPTIONS  IN  EFFECT* 
•OPTIONS  IN  EFFECT • 


NAME  = MA IN, OPT  = 01 ,L1NECNT  = 60. S./E  = u00 .K . 

SOURCE  .EHCUIC.NOL  1ST  .NOOECK  ,LOAO.NLM<’r  ,r.  r o ’ T . 1 0 . NOXREF 


•STATISTICS*  SOUWCE  STATEMENTS  = 27  .PWOORAM  SIZE  = 135* 

•STATISTICS*  NO  DI AGNOST ICS  GENERATED 

• •••••  End  of  COMPILATION  •*•*•*  7 IK  BYTES  OF  CORE  NOT  USEO 


BEST  AVAILABLE  COPY 


37 


6.6 


INTQ  SUBROUTINE  DESCRIPTION 


FUNCTION:  Integrates  both  the  transition  matrix  and  its  inverse  from 
T to  T+H.  [The  subroutine  for  evaluating  the  differential 
equations  is  passed  in  via  the  calling  list.] 

INPUT: 

PS  6x4  array  containing  the  four  intermediate  values  of  the  state 
vector  needed  for  the  Runge-Kutta  algorithm 

PHI  6x6  array  containing  transition  matrix  or  inverse  at  time  T 
T Time  at  beginning  of  interval  of  integration 
H Integration  step  size 
X 9 dimensional  state  vector  at  time  T+H 

FQ  Procedure  which  is  called  to  evaluate  the  right  hand  side  of  the 
differential  equations 

OUTPUT: 

PHI  6x6  array  containing  the  transition  matrix  of  its  Inverse, 
dependent  upon  the  way  subroutine  is  called,  at  T+H 


IL 


LEVEL  21.7  ( JAN  73  ) 0S/360  FORTRAN  M 

COMPILER  OPTIONS  - NAMt=  M4IN,OPT=01 .LlNtCNT=60.SIZE=O0n0K. 

SOURCE .EbCO I C .NOL 1ST .NODECK, LOAD. NO«AP,NOEOIT. ID. NOXREF 
1SN  0002  SUtt  R Oil  T IF'E  INTO(PS.pkI.PmID.T.H.A.FU) 

c suhrout i me  into  integrates  either  transition  matrix  from  t to  t»n. 
C INF  computation  of  The  DERIVATIVES  is  supplied  by  an  external 
C PROCEDURE  AND  IS  DESIGNATED  IN  THE  PARAMETER  LIST  AS  FO.  FS  AWE  THE 
C FOUR  VALUES  OF  X NttDEO  FOP  THE  RUNoE  KUTTA  INTEGRATION  ALGORITHM.  THE 
C VALUE  OF  X AT  T Hi  INTERVAL  OF  INTLOPATION  IS  DESIGNATED  AS  X.  IT  IS 
C USED  TO  EVALUATE  PhIO  AF TER  Tut  VALUt  OF  Pnl  »T  THE  END  OF  TmE 
C INTERVAI  HAS  bEEN  COMPUTED. 


ISN 

0003 

implicit  real*8< a-h.o-zi 

ISN 

0009 

DIMENSION  PS (6.9)  .PHI (6.6) .PH ID (6.61 »S <6.61  .0(6,6) .0 (6*6) .V (6) . 
K C ( 3 ) ,0(9) ,X<1) 

ISN 

0005 

DATA  C/2«0. SJ0.1D0/, 0/100, 2«2D0, 100/ 

ISN 

0006 

00  5 1=1.6 

ISN 

0007 

00  5 J=1 .6 

ISN 

0006 

U< I , J) -PH  I ( I, J) 

ISN 

0009 

5 

S ( I . J 1 =PH I ( I , J) 

ISN 

0010 

DO  30  1=1.9 

ISN 

0011 

IF(I.EO.l)  GO  TO  15 

ISN 

0013 

K*I-1 

ISN 

0019 

DO  10  J= 1,6 

ISN 

0015 

00  10  L = 1 .6 

ISN 

0016 

in 

U(J.L)=S(J.L>*C(K)«Q<J,L) 

ISN 

0017 

15 

CONTINUE 

ISN 

0018 

00  16  J= 1 , 6 

ISN 

0019 

16 

Y ( J) =PS ( J. I 1 

ISN 

0020 

CALL  FmatRX(Y) 

ISN 

0021 

CALL  FO(U.O) 

ISN 

0022 

DO  17  J= 1 , 6 

ISN 

0023 

DO  17  L = 1 .6 

ISN 

002A 

1 7 

0 ( J.L 1 =H*Q ( J ,L ) 

ISN 

0025 

DO  25  J=1  .6 

ISN 

0026 

DO  25  L = 1 .6 

ISN 

0027 

25 

RHl  ( J.L)  =PHl  ( J.L) *D( I ) «0( J.L1/6D0 

ISN 

0026 

30 

continue 

ISN 

0029 

CALL  FMATRX(X) 

ISN 

00  30 

CALL  FOIPHI.PHID) 

ISN 

0031 

return 

isn 

00  32 

eno 

•OPTIONS  IN  EFFECT*  NAME*  «A I N .OPT = 0 1 .L I NE CNT =60 . S I ZE=0000K , 

• OPTIONS  IN  EFFECT • SOURCE .EHCDIC.NOLlST , NOOECk .LOAD , NOMAP .Nr  EDIT. ID. NO X REF 

•STATISTICS*  SOURCE  STATEMENTS  = 31  .PROGRAM  SI7E  = 2008 

•STATISTICS*  NO  DIAGNOSTICS  GEnERmTEO 


END  OF  COMPILATION 


71k  bytes  of  cope  not  used 


OATE 
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6.7 


DPHINV  SUBROUTINE  DESCRIPTION 


EUNCTION:  Computes  derivative  of  the  Inverse  of  the  transition  matrix 
INPUT: 

Q 6x6  array  containing  Inverse  of  transition  matrix 
/FHIM/ 

F 6x6  array  containing  coefficients  of  differential  equations  tor 
inverse  of  transition  matrix 

OUTPUT : 

DQ  6x6  array  containing  derivatives  of  the  inverse  of  the  t.an-.itit 
matrix 

ENTRY  DPHI 


FUNCTION:  Computes  derivatives  of  transition  matrix 
INPUT: 

Q 6x6  array  containing  transition  matrix 
/PHIM/ 

F 6x6  arra/  containing  coefficients  of  differential  equations  for 
transition  matrix 

OUTPUT. 

DQ  6x6  array  containing  derivatives  of  transition  matrix 


DPHINV  FLOWCHART 


COMPUTE 

DERIVATIVE  OF  INVERSE 
OF  TRANSITION  MATRIX 


ENTRY  D P H I 


COMPUTE 
DERIVATIVE  OF 
TRANSITION  MATRIX 


LEVEL 

/l. 7 

( JAN  73  ) 

COMPILER  OPTIONS  - NAMt  = ma[N,OPT  = 0 
SOURCE. EBCDIC. NOL 

1 SN 

0002 

SUBROUTINE  UPH I N V ( 0 . DO ) 

ISN 

0003 

Implicit  real»6 ia-h, o-z> 

! SN 

GOO*. 

CPu"ON/Pm I m/E (6.6) 

ISN 

0005 

DIMENSION  0(6.61.00(6.6) 

ISN 

0006 

DO  10  1=1,6 

ISN 

0007 

DO  10  J= 1 . 6 

ISN 

0006 

DO ( 1 , J) =000 

ISN 

0009 

DO  10  K=  1,6 

ISN 

0010 

10  D0(1.J»=D0(1.J)-0(I.K),F( 

ISN 

0011 

RETURN 

ISN 

0012 

ENTRY  DPHI (O.DO) 

ISN 

0013 

DO  20  1=1,6 

ISN 

0014 

00  20  J=1 ,6 

ISN 

0015 

DO ( 1 . J) =000 

ISN 

0016 

DO  20  K = 1 , 6 

ISN 

vO  1 7 

20  DO ( I « J) =DQ ( I , J) *F ( I ,K ) »Q ( 

ISN 

0016 

RETURN 

ISN 

0019 

ENO 

OS/360  FORTRAN  H 


•OPTIONS  I N EFFECT • NAM£=  MA I N , OP T =0 1 . L I Nt CN T = 60 . S I ZE =0000K . 

•OPTIONS  IN  EFFECT*  SOURCE ,E6CD I C . NOL I ST . NODlC* . LOAD .NOMAP . N^f 0 1 T , I D . NOXREF 

•STATISTICS*  SOURCE  STATEMENTS  = 16  .PROGRAM  SIZE  * HOO 

•STATISTICS*  no  diagnostics  generated 

• •••••  E NO  OF  COMPILATION  **•*•*  79*  BYTES  OF  COPE  NOT  USED 
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6.3 

FUNCTION: 

INPUT: 

X 

OUTPUT: 


FMATRX  SUBROUTINE  DESCRIPTION 

Computes  matrix  which  relates  transition  matrix  and  Its  inverse 
to  their  derivatives 


9 dimensional  array  containing  position,  velocity  and  acceleration 


6<6  dimensional  array  containing  matrix  which  relates  transition 
matrix  and  its  inverse  to  their  derivatives. 


SUBROUTINES  USED:  GRVPAR,  ROTPAR,  DRGPAR 


FNATRX  flowchart 
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bestavailad::  copy 


If.  VEL  21.7  < JAN  73  ) 


OS/360  FORTRAN  M 


UATE 


COMPILE®  OPTIONS  - NAME*  MAIN.OPT=O1.lINECNT*60,SIZE«0000K. 

SOURCE  *EBCO I C »NOL 1ST .NODECK  .LOAD .NOmAP.NOEOI T • I D.NOXREF 
ISN  0002  SUHROUT I NE  F«ATRX(X> 

C THIS  SUBROUTINE  COMPUTES  THE  F MATRIX  USEU  IN  T-E  OIEFEPENTIAL 
C EQUATION  FOP  PHI  AnD  PHI  INVERSE.  THE  DERIVATIVE  OF  PHI  IS  EQUAL  TO 
C F TIMES  PHI.  THE  Otx I VAT  I VE  OF  PHI  INVERSE  IS  f ''UAL  TO  MINUS  PHI 
C INVFRSE  TIME  PH  1 • THESE  EQUATIONS  ARE  INTEGRATED  TO  OBTAIN  PHI.  F IS 
C STORED  IN  THE  VARIABLE  F. 

ISN  0003  IMPLICIT  REAL*B(A-h,0-2) 

ISN  000*  COMMON/PHIM/F (0,61 

ISN  OOOS  DIMENSION  X<9) 

ISN  0006  DO  5 1=1.3 

ISN  0007  DO  S J=  1 * 6 

ISN  0008  F ( I , J) =000 

ISN  0009  IFU.EO.IO)  F(I,J)=1D0 

ISN  0011  S CONTINUE 

ISN  0012  DO  10  I =*.6 

ISN  0013  DO  10  J=  1 .6 

isn  ooi*  in  f ( i . j>  =ooo 

ISN  001S  CALL  GRVPAR(X) 

ISN  0016  CALL  POTPAR 

ISN  0017  CALL  ORGPAR ( X 1 

ISN  0018  RETURN 

ISN  0019  ENO 


•OPTIONS  IN  EFFECT* 

•OPTIONS  |N  EFFECT* 

•STATISTICS*  SOURCE  STATEMENTS  * 
•STATISTICS*  NO  DIAGNOSTICS  GENERATED 


NAME  = ma I N .OPT  = 01 .L InECNT  =60.SIZE=OOOOK. 

SOURCE .E0CUIC.NOL IST.NCPtCr « LOAD . NOMAP . NOE D I T . ID .NOXREF 
16  .PROGRAM  SIZE  ■ SOB 


t NO  OF  COMPILATION 


7VK  BYTES  OF  CORE  NOT  USEO 


6.9  DRGPAR  SUBROUTINE  DESCRIPTION 

In  this  application  this  subroutine  is  a dummy  subroutine.  Drag  is  not 
modeled  in  either  the  powered  flight  or  free  fall  segments  of  the  trajectory. 
However,  when  the  reentry  portion  of  the  trajectory  is  considered,  these 
parti als  will  have  to  be  computed. 
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■sm — wm 


-wwiP^r 


BESI  AVAI'J"'--  COPY 


LEVEL  Z1.7  l JAN  73  ) 

CC-PlcER  OPTIONS 


ISN 

r0b<! 

C 

C 

c 

SUBROUTINE  IIOGPAR(X) 
SU-POUTI'.'f  LRGPaR  COMPUTES  t 

*ITh  respect  to  position  ano 
F. 

I SN 

oOO  3 

Implicit  reac'hia-n.o-Z) 

! SN 

■ 3u  . 

COM«.  ON/pm t m/F  16.61 

I SN 

GoOs 

DIMENSION  .(V) 

ISN 

0006 

S=1D0 

ISN 

0007 

RE  TURN 

I SN 

uOOo 

END 

OS/3bO  FORTRAN  H 

M‘**c  = MAlM.OPT^r  1 ,L  lNECNT=bO.SI  Zt  =0000K  , 

SOURCE  .EnCU  IC  .NOL  1 S f , NOUECk  ,cOAu  . NO*.'  AP  , NOEO 1 T , IO.nOaREF 


"•OPTIONS  IN  EFFECT' 
•OPTIONS  IN  FFFECT' 


•STATISTICS*  SOURCE  STATEMENTS  = 

•STA* 1ST  ICS*  NO  DIAGNOSTICS  generated 


name ; MAIN«0PT-01.LINECnT-60.SIZE-0000K» 

SOUPCE .EbCulC .NOL 1ST .NOUlCK .LOAD . NOMAP . N',E u I T . ID.NOXREF 
7 .PROGRAM  SIZE  = >?38 


UATE 


End  of  CGMPICATION 


7Rk  bYTES  OF  CORE  NOT  USt'O 


6.10 


ROTPAR  SUBROUTINE  DESCRIPTION 


FUNCTION:  Computes  partial s of  acceleration  due  to  the  rotation  of  the 
coordinate  system  with  respect  to  position  and  velocity 

OUTPUT: 

F Partials  of  rotational  acceleration  with  respect  to  position 
and  velocity  are  added  to  6x6  dimensional  array  containing  F 

SUBROUTINES  USED:  None 


BtSl  AVM IABI£  COP>f 


Lf  VEL 

<n ./ 

I SN 

GOOZ 

1SN 

000J 

ISN 

000a 

ISN 

COOS 

ISN 

0006 

I SN 

000/ 

ISN 

0008 

ISN 

ooov 

ISN 

ooio 

ISN 

0011 

ISN 

ooia 

OS/ 360  F0PT3AN  M 


UATt 


LOmPHEW  OPTIONS 


NAME  - MA  IN,  O^T  = 01  • LINECNT=60«SIZE  = /'000k. 

SOUwCE.LtiOOiC.NOL  1ST  .NODE  O'  . LOAU  . NO'' AP  . NOLU  I T, IU.NOXHEF 
SUPOOUTINF  OOTPAW 

C SUBPOUTInE  POT  Paw  COMPUTES  T Ht  PARTIAIS  Of  AlCFLEBaT  I ON  DUE  TO  TmE 
C ROTATION  Of  The  EARTH  with  PESPlCT  TO  POSITION  • Nu  VELOCITY.  THESE  ABtC 
C AOOEO  INTO  THE  f MATRIX. 

IMPLICIT  WEAL*«<A-H,0-Z> 

COMmqN/PhIm/F (6.6) 

OAIA  W/0. 7^S1 IS1D-A/ 
w2=w*w 

F (A. 1 ) =F (A, 1 I 
F (S.?>  =F(S.2) *W a 
f <A,S>  =F  <A,S> ♦w*ZDO 
F (S. a)  =F  (5. A)  -v**00 
RETURN 
* END 


•OPTIONS  IN  EFFECT* 
•OPTIONS  IN  EFFECT* 


•STATISTICS*  SOURCt  STATEMENTS  * 
•STATISTICS*  no  diagnostics  generated 

••••••  END  OF  COMPILATION  •••••• 


name  = MA IN. OPT  = 01 .LINECNT=60.SIZE=OOOOK , 

SOURCE  .EFiCOIC.NOLIST.NODtCK  .LOAD.NOMAI'.NOEOI  T.  ID.NOXNEF 
11  .PROGRAM  SIZE  = </G0 


7VK  BYTES  OF  COPE  NOT  USED 
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6.11 


GRVPAR  SUBROUTINE  DESCRIPTION 


FUNCTION:  Computes  partial s of  gravity  acceleration  with  respect  to  position 
INPUT: 

X 9 dimensional  state  vector.  Only  position  terms  are  used 

/TIG/ 

IN  Flag  to  indicate  whether  the  partial s are  to  be  computed 
numerically  or  analytically 

OUTPUT: 

F 6x6  array  containing  gravity  partials 
SUBROUTINES  USED:  EGRAV , FVPMAT 

COMMENIS:  GRVPAR  must  be  called  before  ROTPAR  or  DRGPAR.  IN  = 0 gives 
numerical  partials.  The  analytic  partials  and  the  gravity  is 
based  upon  the  J and  D model. 
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BEST  AVAI 


i r 

ui  IU  — 


v * t r 


( JAN  / 3 t 


OS/J60  FOOT-TAN  M 


LO-PILl*  OPTIONS  - NA**t  = MAIN.OPT  = 01.LlNECNTa60.SIZE=OOOOK, 


1SN  t.  0t>  «• 


sr.-p  CE  . e -3C  U ! C • NUL  i 'll  .NODFCk  .LOAU.NO-AP.NOCOI  T , IO.nOaPEF 


Sii-hOOT  1 “ 

S-»--r.»  - T 1 nE  G. 
CjJ«v1T»  * I Tn 


> I.WVPAPIX) 

ui’*-'  COMPUTES  TMf.  ► 
i-tSPECT  TO  POSIT  10'. 


* J T 1 AL  c OF  ACLtl  f-»T  ION  UUt  TO 
AN0  VELOCITY.  TmESE  Ahl  ADDED  INTO 


T Hi 


•A  1-1  * F , 


ISN 

oOL'3 

implicit  ueal»s(a-m,o-Z) 

T <;s 

1 OCA 

LOmmON/Ph I M/F (fc.b) 

I <;►' 

v.  uos 

COM'")N  /til/  In 

1 

1 0G6 

Tm-  * SION  A*  ( J)  ,Y  (6)  ,AY  ( 

l C*N 

0 00/ 

. - 1 «E  M'  ION  G(  J.3) 

! s* 

COOh 

IE  I Iv.E'l.O*  SO  TO  TO 

! ^N 

0 10 

On  10  1-1.3 

* SN 

00  1 1 

l-i)  10  1=1.6 

1 s»' 

"OIP 

v ( J)  **  ( 3) 

• c » 

'01  * 

10 

Y(J1 =*(J1 

» <;* 

OlA 

< 1 1 > =y ( 1) .SOO 

t r • • 

01S 

V ( I 1 =V  ( 1 1 -SIlO 

1 r-.N 

01E 

cam  * r,w a v t y . a y> 

! r*" 

101/ 

CALI  E OP  A V ( V . A V ) 

I S»J 

OOln 

DO  PO  J= 1 , 3 

ISM 

O01T 

K = j.  3 

I VJ 

00P0 

PG 

f (a  ,i | = (AY(J|-AV(J;  1/101 

* *N*. 

(0*1 

30 

CONI INUE 

! f>N 

LOPP 

m 

( ON  T I NUF. 

I SN 

Ou«r  J 

1*  ( IN.f  0. 1 ) GO  TO  60 

I SN 

OOPS 

(ALL  FvpMAflA.G) 

I SN 

o 0/6 

DO  SO  !=l,J 

I SN 

00  c 1 

no  So  J-  1 • 3 

I SN 

OOPH 

* = I ♦ 3 

I SN 

00/* 

so 

F<«. J»S Gil fj) 

I SN 

00  JO 

to 

C(#AlJNUE 

I SN 

00  31 

wtV'0* 

: v 

C 0 JP 

E NO 

N 

T f r*NS  |n 

E F E t C T • 

i|AME  r MA  | N.OP  f = 0 1 tl  1 

f I*r 

s IN 

E‘FfCT» 

StlUPCf  .EnCNIC. NOLIST, 

■r-' 


•STATISTICS"  SOURCE  STATEMENTS  = 

•STATISTICS*  NO  0 I At-NOS  TICS  OE  NEGATE  0 


IN*  <.N  I - fsO.Sl  /t s OOOOf  . 


11  .PPOUPAM  S 1 7E  * >/3<* 


End  uf  compilation 


V)K  OYTtS  OF  COPt  NOT  USED 


DATE 
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A 


6.12  EGRAV  SUBROUTINE  DESCRIPTION 

FUNCTION:  Computes  gravity  acceleration  using  J and  D model  of  gravity 
INPUT: 

Y 9 dimensional  array  containing  state  vector.  Only  first  3 
elements  are  used 

OUTPUT: 

GY  3 dimensional  array  containing  acceleration  for  three  variables 
SUBROUTINES  USED:  None 

COMMENTS:  For  descripton  of  gravity  mathematics  see  [4]. 
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A 


LEVft  2 1.7  « JAM  13  1 


UtTI 


BEST  AVAILABLE  COPY 


OS/360  FOftTRAN  M 

COMP  I L t ° OPTIONS  - N»Mt  = M A I N .OPT  =0  1 .L  I Nt  CNT  = 60  • S I It  = 0000*  • 

SOUPCE «E  HCl)  1 C «NOL 1ST .NOUtCK.COAO.Ntl  (AP.NOtOI  T.  1D.N0K0EF 
ISM  0002  SiiPPOUTIMf  tGPAV(Y.GY) 

ISM  OOOJ  IMPLICIT  PtAL«d  (A-H.O-ZI 

ISM  000*.  U I MENS  I ON  Y(l).GYU) 

ISM  000*5  OATA  GM,  CllNO  .COM  A ,C  JC  AS/ 1 .9076468S  016. 

• 0.1027VS0G  0-4,  20926672.6  DO.  0.7H1SS2I9  01?/ 

c»*»«»*»**»« ••••••••*•••••••••••••••••••••«••••••••••••••*•••«• ••••••••• 

C COMPUTt  THE  ACCtLEPAT 1 ON  IN  EE 6 Out  TO  GRAVITY 

C»  This  subroutine  IS  TAKEN  FPOm  A TRW  REPORT  on  JOT  SIMULATION  • 


ISN 

(1006 

A * COMA 

\t  « <*  « 1 

ISM 

0007 

0 * COMO 

ISM 

0006 

P?  * 0. 

ISN 

0009 

DO  10  1*1.3 

ISN 

0010 

10  *2  * R?  » Y(I)*Y(I) 

ISN 

0011 

W4  = P?»«? 

ISN 

001? 

/2P2  = Y(31»Y(3»/P? 

ISN 

0013 

A?  * A • A 

ISM 

0019 

DAP  = D*A2«A2 

ISM 

001S 

h s (GM/P?I  • ( 1 .00  ♦ (CJCAS/K2) * ( 1 .00  - S.D0*7?(<?»  ♦ <0A4/P«.;» 

ISM 

0016 

* ( 0 . 4266 7 14  300  - Z2P2M6.U0  - 9 .00»Z2M2H  » 

A?  = OSORT <R2) 

ISM 

0017 

Gy ( 3>  = -(Y(J)/A?)»(H  . (0M/P9) •<?.00*CJCAS  ♦ (Daw/P?) • 

ISN 

0016 

• <1.  71428671400  - 9.D0*Z?P?m 

GY ( 1 ) = -<H»Y<1>/A2) 

ISM 

0019 

GY  1 2 > = -1H*Y12>/A2) 

ISM 

0020 

RETURN 

ISN 

00c  1 

t MU 

•OPTIONS  IN  EFFECT*  NAME  * MA IN .OPT=0 1 .L lNECMT=o0 .SI Zt*OOOOK , 

•OPTIONS  IN  EFFECT*  SOUPCE. EBCDIC. NOLIST. NOOtCK.LOAO. nomap, M^EOIT. ID. NOXBEF 

•STATISTICS*  SOUPCE  STATEMENTS  = 20  . PROGPAM  SIZE  * 814 

•STATISTICS*  no  DIAGNOSTICS  generated 

• •••••  End  UF  COMPILATION  ••••••  79k  bytes  of  COPE  MOT  USE 0 


6.13 


DIFFUN  SUBROUTINE  DESCRIPTION 


FUNCTION:  Computes  acceleration  due  to  gravity  and  acceleration  sensed 
because  of  rotating  coordinate.  Computes  velocity  of  body  in 
motion 

INPUT: 

N Number  of  differential  equations  (never  used  because  It  is 
always  6) 

T Time  (not  used) 

Y 9 dimention  array  containing  position,  velocity  and  acceleration. 
Only  position  and  velocity  are  used 

OUTPUT: 

DY  6 dimensional  array  containing  gravity  and  rotation  accelerations 
and  velocities. 

SUBROUTINES  USED:  EGRAV 


1 


1 


BESrAVAILACi  COPY 


LEVEL  21. 7 C JAM  73  ) 

COMPILED  OPTIONS 

ISN  0002 


0S/360  FORTRAN  M 


DATE 


ISN  0003 
ISN  O00* 


ISN  0005 
ISN  000b 
ISN  0007 

isn  oooa 

ISN  0009 
ISN  0010 
ISN  0011 
ISN  0012 
ISN  0013 
ISN  001* 
ISN  00  lb 
ISN  001b 


NA*t=  MA1N.OPT«01 .LlNECNT=b0.SIZt«0000K. 

SOURCE .EbCOIC.NOLlSI . NODE CK, LOAD. NOvAP, NOEOIT. 10. NOXREF 
SUBROUTINE  DIFFUNIN.T.Y.DY) 

C COMPUTE  THE  ACCELERATION  Out  TO  GRAVITY  IN  EFG  COORDINATES 
IMPLICIT  PtAL»b  (A-H.O-Z) 

DIMENSION  Y(l> .OYC1) .GYI3) 

(;•••••••••••• •••••••••••••••••••••••••••••••••••••••••••••• 

C»  m IS  EARTH'S  ANGULAR  VELOCITY  IN  PAOI AnS/SECONO  • 

C*«« •••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• 

v = 0. 72921 1510-* 

VSO  = R«W 
*>?  = ?.00»R 
20  OY(l)  s VIA) 

OY (2)  = V(5) 

OY 13)  s y I b ) 

CALL  FGRAVIY.GY) 

OYI4)  = GYll)  ♦ *2*Y<5)  * WSO*YU) 

OY 15)  = Gr (2)  - *2»Y(A)  ♦ *SQ«Y<2) 

DYIb)  * GT ( 3) 

RETURN 

ENG 


•OPTIONS  IN  effect* 
•OPTIONS  IN  EFFECT* 


•STATISTICS*  SOURCE  STATEMENTS  * 
•STATISTICS*  NO  OIA&nOSTICS  GENERATED 

••••••  END  OF  COMPILATION  •••••• 


NAME  s MAIN.OPT  = 01 .LlNECNT=bO.SIZE*OOOOK, 

SOURCE .EbCOIC .NOL I ST .NOOECA . LOAO. NOMAP. NoEDI T . ID.NOXREF 
15  .PROGRAM  SIZE  « *96 


79X  BYTES  OF  CORE  NOT  USEO 


6.14 


FVPMAT  SUBROUTINE  DESCRIPTION 


KUNCTJ_ON:  Computes  partial s of  gravity  acceleration  with  respect  to  position 
INPUT: 

X 9 dimensional  array  containing  position,  velocity  and  acceleration. 
Only  position  is  used 

OUTPUT : 

F 3x3  array  containing  partials  of  gravity  acceleration  with  respect 
to  position 

SUBROUTINES  USED:  None 

COMMENTS:  Partials  are  computed  for  J and  0 model  of  gravity 
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COMPUTE 


PARTIALS  OF  ACCELERATION 
WITH  RESPECT  TO  POSITION 


OA  »€ 


BEST  AVAILABLE  COPY 


IIVFL  21.7  ( JAN  7 J > 


OS/J60  FOPT9AN  H 


ISN 


ISO 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

I«N 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 

ISN 


COMPILER  OPTIONS  - NAMt = ma;n,opT=01 .1 INtCNT=60.SIZE=0000K. 

SOUMCt .E9CDIC.N0L I Sf .NOUECK .LOAO.NOmap .NOE  0 1 T , 10 .NOXPEF 
0002  SUBROUTINE  EVPMAT(X.F) 

C* 

C*  COMPUTE  THE  PAWTIALS  OF  ACCtLEBAT ION  x 1 TM  PESPECT  TO  POSITION 

C* 

C*  PI  = P(  p.  (l./P*»2),  <1./M**9>  ) / P(  X,  Y.  Z) 

C* 

c*  p?  = p(  i / p<  x.  v « z ) 

c* 

C*  P3  = P(  x/M,  Y/P,  Z/R)  / P<  X,  Yt  Z I 

C* 

C*  P9=P(H)/P(X.Y.  Z> 

C* 


0003 

000** 

0005 

0006 
0007 
0006 

0009 

0010 
0011 
0012 
0013 

0019 

0015 

0016 

0017 

0018 

0019 

0020 
1/021 
0022 
0023 
0029 

0025 

0026 
0027 
0026 

0029 

0030 

0031 

0032 

0033 
0039 

0035 

0036 

0037 

0038 


0039 

0090 


IMPLICIT  REALMS  (A-h.O-Z) 

PEAl«8  OM.JA2 

DIMENSION  X( l ) .PI (3,3) .P2( 31 .P3(3»3) «P9(3> ,F(3.3> 

DATA  CONO/  0.10279500  0-9  / 

DATA  CJCAS/  0.711155219  012  ✓ 

UATA  ¥/  0.72911510-9  / 

DATA  GM/  1 .90769665  016  / 

DATA  A/  20925672.6  OO  / 

EQUIVALENCE  (CJCAS. JA2) , (0. CONO) 

P2  = X ( 1 ) *X ( 1 ) ♦ X (2) *X (2)  ♦ X (3) *X (3) 

» = OSOPT (B2> 

P3  a B*P2 
P9  = R2*B2 
96  = P3*P3 
UO  10  1=1.3 

Pl(l.I)  = X(I> /R  . 

PI (2.1)  = -2.U0*X(l)/R9 
10  PI (3.1)  = -9.00*X(I)/B6 
Z2  = X ( 3) ®X ( 3) 

P2  ( 1 ) = -2.U0*X(1)»Z2/P9 
m2 ( 2)  = -?.u0*x(2) *72/99 
p? < 3)  = 2.00*x<Ji*(»2  - Z2I/P9 
P3I1.1)  = (P2  - X ( 1 ) *x ( 1 ) )/P3 
P3I1.2)  = -x ( 1 ) «x (2) /R3 
93(1.3)  = -X(1)»X(3)/R3 
P3I2.1)  = P3I1.2) 

P3 (2.2)  = (P 2 - x (2>*x (2) )/H3 
93(2.3)  = -x (2)*X(3)/»3 
93(3.1)  = P3I1.3) 

93(3.21  = P3 ( 2 .3) 

P3 ( 3. 3)  = (P2  - Z2I/P3 
72  = Z2/R2 
0 A9  = 0*A*A*A*A 

H = (Gm/P2)*(1.OG»(JA2/p2)*i1.D0  - 5.uO*Z2>  • (UA9/B9)* 

1 (0.9265719300  - 72* (6.00  - 9.00*Z2>>) 

00  20  1=1.3 

P9 ( I ) = <GM/H2)*(  ( JA2/P2) • ( -5.00*P2 ( 1 1 ) ♦ JA2* ( I .00  - S.D0*Z2) 

1 *91(2.1)  - (DA9/«9)*(Z2*(-9.O0*P2(I) ) * (6.00  - 9.00*72) 

2 *P2 ( I ) ) ♦ OA9*(0. 92 o5 719300  - Z2*(6.U0  - 9.00*Z2)) 

3 *P1 (3. I ) ) ♦ M*R2*P1 (2,  II 
20  CONTINUE 

GZ  = -(X(3)22)*(H  ♦ (GM/P9! *(2.O0*JA2  ♦ (0A9/P2) 
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BESI  AVAICABLE  COPY 


•U. 71428571400  - 4.U0»Z2)>» 


ISN 

0041 

DO  30  1-1.3 

ISN 

004* 

F(l.I)  = -H»R3 ( 1 • I ) - (X(1)/H)«P4(I) 

ISN 

0043 

F(2.I)  = -M*°3 (2.1)  - ( A ( 2 I /R ) *P4 ( I ) 

ISN 

g044 

F ( 3 » I ) = - (X (3) /R) * (P4 ( I ) . ( GM/R4 ) * 

1 ♦ UA4* ( 1 . 7 1 4285 7 1400  - 4.00 

2 ♦ Gh*(2.00»JA 2 ♦ (0A4/W2) * ( 

3 •PII3.Il)  » (R«GZ/X(3) )»P3C 

ISN 

004S 

30  CONTINUE 

ISN 

0046 

•SO  * 

ISN 

0047 

RETURN 

ISN 

0048 

ENO 

*.D0«Z2)  ) 


•OPTIONS  IN  EFFECT*  NAME  = MA I N ,OPT=0 1 .L I NE CNT =60 . S I 7E*0000K , 

•OPTIONS  IN  EFFECT*  SOURCE .EflCO I C . NOL I ST , NODECK . LOAD .NOMAP , NOEUI T . ID.NOXREF 

•STATISTICS*  SOURCE  STATEMENTS  « 47  .PROGRAM  SIZE  * 1BS4 


•STATISTICS*  NO  OIAGNOSTICS  GENERATED 

••••*•  EMU  OF  COMPILATION  *•••*•  67K  BYTES  OF  COPE  NOT  USED 

•STATISTICS*  NO  DIAGNOSTICS  THIS  STEP 
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6.15 


CON  SUBROUTINE.  DESCRIPTION 


! 


FUNCTION: 


INPUT: 

/VAR/ 

H 

XO 

X 

PHIO 

PHI 

PHIDO 

PHID 

PH  I NO 

PHINDO 

PHIN 

PHIND 

ID 

OUTPUT: 

/VAR/ 

XC 


PHIC 


Computes  constants  for  the  Spline  interpolating  polynomials  for 
the  position  in  the  state  vector  and  for  the  transition  matrix 
and  its  inverse 


Step  size  of  integration 

9 dimensional  array  containing  position  and  velocity  at  the 
beginning  of  the  interval 

9 dimensional  array  containing  position  and  velocity  at  the  end 
of  the  interval 

6x6  dimensional  array  containing  transition  matrix  at  the 
beginning  of  the  interval 

6x6  dimensional  array  containing  transition  matrix  at  the  end 
of  the  Interval 

6x6  dimensional  array  containing  derivatives  of  PHIO 

6x6  array  containing  derivatives  of  PHI 

6x6  array  containing  the  inverse  of  PHIO 

6x6  array  containing  derivatives  of  PHINO 

6x6  array  containing  the  derivatives  of  PHIN 

6x6  array  containing  the  derivatives  of  PHIN 

FLAG  TO  COMPUTE  COEFFICIENTS  FOR  SPLINES  FOR  TRANSITION  INVrRSE 


6x3  array  containing  coefficients  for  the  qulntlc  Splines  for 
position.  The  first  Index  indicates  the  polynomial  coefficient 
and  the  second  indicates  the  variable 

4x6x6  array  containing  coefficients  for  the  cubic  Splines  for 
the  transition  matrix.  The  first  Index  indicates  the  coefficient 
and  the  second  and  third  indices  indicates  the  matrix  element. 
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PHINC  4x6x6  array  containing  coefficients  for  the  cubic  Splines  for  the 
inverse  of  the  Transition  matrix.  The  Indcles  are  the  same  for 
RHIC 
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